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Abstract 

Hypoxia in Chesapeake Bay has substantially increased in recent decades, with detrimental effects on macrobenthic 
production; the production of these fauna link energy transfer from primary consumers to epibenthic and demersal 
predators. As such, the development of accurate predictive models that determine the impact of hypoxia on macrobenthic 
production is important. A continuous-time, biomass-based model was developed for the lower Rappahannock River, a Bay 
tributary prone to seasonal hypoxia. Phytoplankton, zooplankton, and macrobenthic state variables were modeled, with a 
focus on quantitatively constraining the effect of hypoxia on macrobenthic biomass. This was accomplished through 
regression with Z': a sigmoidal function between macrobenthic biomass and dissolved oxygen concentration, derived using 
macrobenthic data collected from the Rappahannock River during the summers of 2007 and 2008, and applied to compute 
hypoxia-induced mortality as a rate process. The model was verified using independent monitoring data collected by the 
Chesapeake Bay Program. Simulations showed that macrobenthic biomass was strongly linked to dissolved oxygen 
concentrations, with fluctuations in biomass related to the duration and severity of hypoxia. Our model demonstrated that 
hypoxia negatively affected macrobenthic biomass, as longer durations of hypoxia and greater hypoxic severity resulted in 
an increasing loss in biomass. This exercise represents an important contribution to modeling anthropogenically impacted 
coastal ecosystems, by providing an empirically constrained relationship between hypoxia and macrobenthic biomass, and 
applying that empirical relationship in a mechanistic model to quantify the effect of the severity, duration, and frequency of 
hypoxia on benthic biomass dynamics. 
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Introduction 

Macrobenthic organisms (retained on a sieve size >500 urn) are 
of importance to ecological processes in estuarine ecosystems like 
Chesapeake Bay [1], regulating or modifying most physical, 
biological, chemical, and geological processes [2]. Macrobenthos 
influence these sediment geochemical and physical properties [3] 
through bioturbation, the biological reworking of sediments [4] . In 
the estuarine environment, macrobenthos are a primary pathway 
through which organic carbon is cycled in the sediments [1]. 
However, the role of macrobenthos is limited to normoxic 
conditions; anaerobic metabolism becomes an increasingly 
important means of recycling organic carbon during hypoxia 
[5]. Macrobenthos also serve as a major energetic link between 
primary producers and demersal fish and epibenthic predators [6] ; 
however, the sessile nature of macrobenthos makes them 
susceptible to natural and anthropogenic perturbations such as 
hypoxia [7], a significant concern given the documented 
importance of estuarine macrobenthic communities [1]. 



Since colonial times, the number of humans in Chesapeake Bay 
watershed has grown exponentially, with a 3-fold increase during 
the last 100 years [8]. Human activity adversely affects land 
topography, chemistry of the Earth's atmosphere and water, rates 
and balance of biogeochemical processes, and biodiversity [9]; 
Chesapeake Bay estuary is no different. Anthropogenic distur- 
bance has greatly increased the flux of nitrogen and phosphorous 
compounds through land clearing, application of fertilizer, 
discharge of human waste, animal production, and combustion 
of fossil fuels, leading to eutrophication of the Bay [10]. Hypoxia, 
dissolved oxygen (DO) concentrations <2 mg 0 2 1 [11], is 
closely associated with eutrophication, an increase in the rate of 
supply of organic matter to a system [12]. Low DO concentrations 
have been documented in mainstem Chesapeake Bay since the 
early 1930s [13] and in the Potomac since the 1910s [14]. 
Presently, seasonal hypoxia forms in early to late spring and lasts 
approximately 120 days, with the most severe low DO events 
occurring during mid-summer in mainstem Chesapeake Bay [8]. 
From the 1950s to the present, hypoxic volume in the Bay has 
increased substantially, from approximately 3 km to 10 km' [15]. 
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This increase is of concern given documentation of low DO 
impairing growth and reproduction and stressing living resources, 
increasing faunal susceptibility to disease and other environmental 
stressors [11,1 6-25] . 

As hypoxia continues to increase in the Bay and many other 
coastal systems worldwide [15,26], the development of accurate 
predictive models that quantify the ecological impacts of hypoxia 
becoming increasingly important. Recent models have begun to 
include functions relating hypoxia to increased mortality and/ or 
reduced filtration or ingestion of macrobenthos, fish, and shellfish 
(e.g. [27-30]), but typically these functions are hypothetical and 
have yet to be constrained by empirical data. Models generally 
take two forms; detailed complex models that attempt to replicate 
as much of the natural environment as possible [28,31], and more 
simplistic models that only incorporate what is functionally 
necessary [32]. The focus of the current study was to model the 
effect of observed hypoxia, and various scenarios of hypoxia, on 
the biomass of macrobenthos; therefore our approach was to force 
DO concentrations into a simple model to predict macrobenthic 
responses, rather than to first simulate hypoxia with a highly- 
resolved, linked hydrodynamic -water quality model. The ap- 
proach we developed here can be incorporated into more complex 
models as they continue to add additional mechanistic detail with 
respect to ecological responses to hypoxia. We began by 
empirically constraining a formulation for hypoxia-induced 
mortality, incorporated this function into a simplified version of 
the Chesapeake Bay Eutrophication Model [28], and used the 
resulting model to predict the potential effects of hypoxia on 
macrobenthos in the Rappahannock River, a tributary of 
Chesapeake Bay that experiences seasonal hypoxia. 

Materials and Methods 

2.1 Study Area 

Seasonal hypoxia occurs throughout Chesapeake Bay and some 
of its tributaries during the summer months [8] , but in the lower 
Chesapeake Bay, the Rappahannock River is the only major 
tributary with the hydrography that allows for the development of 
sustained seasonal hypoxia [33]. In the Rappahannock, a 
combination of tidal mixing and advection of undersaturated 
mainstem waters into the tributary controls the seasonal hypoxia, 
which develops in late May and abates in early September [34,35]. 
Macrobenthic taxa in the Rappahannock River exhibit a variety of 
life history traits, such as resiliency and recruitment, which affects 
their response to seasonal hypoxia [1]. These responses allow some 
taxa to rapidly recruit to areas post-hypoxia, while other species 
struggle to re-establish, influencing variation in macrobenthic 
biomass; information on the functional characteristics of the 
dominant macrobenthos taxa in the lower Rappahannock are 
displayed in Table 1 . No permits were required for the described 
study, which complied with all relevant regulations. 

2.2 Field Collection 

Macrobenthos data from a previous study were used to develop 
a formulation for hypoxia-induced mortality in our benthic model 
[36]. During the summers of 2007 and 2008, two random sites 
were chosen in each year in the lower Rappahannock for 
continuous monitoring of water quality, and biweekly sampling 
of macrobenthos from May to October. Based on water quality 
data from the Chesapeake Bay Long-Term Benthic Monitoring 
Program (www.baybenthos.versar.com) one normoxic site and one 
site known to experience seasonal hypoxia were chosen in each 
year. Data from a hypoxic and normoxic site in both 2007 and 
2008, four sites total, were used for model construction. 



At each monitored location a Hach DS500X Hydrolab 
datasonde was deployed approximately 0.5 m above the sediment 
surface attached to a small tripod. DO concentration, salinity, 
temperature, and depth were recorded every 20 minutes. 
Approximately every two weeks the datasonde was replaced with 
another Hydrolab datasonde, and a sediment sample was collected 
with a Young grab (440 cm 2 to a depth of 10 cm) for benthic 
community analysis. Grabs were sieved in the field through a 
0.5 mm screen, and organisms and detritus retained on the screen 
transferred into labeled jars, preserved in a 10% formaldehyde 
solution, and stained with Rose Bengal. Samples were processed to 
identify and enumerate each species present as described in Dauer 
and Llanso [37]. Ash-free dry weight (AFDW) biomass was 
measured for each species by drying to a constant weight at 60°C 
and ashing in a muffle furnace at 500°C for four hours. 

2.3 Model Construction 

A continuous-time, biomass-based model was constructed using 
STELLA Modeling and Simulation Software®. The model was 
based on the benthic sub-model in the 2002 Chesapeake Bay 
Eutrophication Model [28,38], and contained three governing 
equations. All formulations and parameter values from the parent 
model [28,38] were used in the current project with the exception 
of the modifications described in sections 2.3, 2.4 and 2.5. 
Phytoplankton biomass was modeled as: 



Wa)*P-PR] 



(1) 



where: 

P = phytoplankton biomass (g C m ) 
G = growth rate of phytoplankton (d ) 
R = respiration rate of phytoplankton (d ') 
Wa = phytoplankton settiing loss (d ) 
PR = predation on phytoplankton (g C m J d '). 
Zooplankton was modeled as the combined biomass of micro- 
and mesozooplankton for simplicity as: 



St 



M=[(Gz-BMz-Mz)*M-PRz] 



(2) 



where:M = zooplankton biomass (g C m~ ) 

Gz = growth rate of zooplankton (d ') 

BMz = basal metabolic rate of zooplankton (d ') 

Mz — hypoxic mortality (d ') 

PRz = predation on zooplankton (g C m - ^ d -1 ). 

Macrobenthos were modeled as the combined biomass of 
deposit and suspension feeders, as: 



St 



B- 



ni2 * 10 6 



* (POC+PM) *k m „i *B 



(3) 



-[(r-ff-m) *B] 



where: 

B = macrobenthos biomass (g C m 2 ) 
a = assimilation efficiency for carbon/ 0 



ingestion rate of 



macrobenthos (g sediment g C 1 biomass 1 d ) 

m 2 = sediment solids concentration (kg 1 )POC = sediment 
particulate organic carbon concentration (g C m -i ) 

PM = phytoplankton and zooplankton biomass (g C m~ 3 )k mnl 
= Michaelis-Menton growth limitation term for carbonr = 
respiration rate of macrobenthos (d~ ) 
P = predation rate (m 2 g C d _1 ) 
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Figure 1. Relationship between DO concentration and macrobenthos biomass. Comparison of macrobenthos biomass and DO 
concentration from the four Rappahannock River sites monitored bi-weekly from May to October. Trendline is a sigmoid curve, where equation 
1.97 



Z' = 



1+e 



: — 3.35 
0.15 



doi:10.1371/journal.pone.0084140.g001 



m = hypoxia mortality rate (d ). 

Phytoplankton and zooplankton groups were included in the 
model given the tight benthic-pelagic (B-P) coupling that exists in 
estuarine and shallow coastal systems [1,39-41], and the 
importance of both groups as a source of food for macrobenthos 
[41]. The recycling of nutrients is another major feature of B-P 
coupling [42], but nutrient recycle terms are not included in the 
model to maintain simplicity. Our model excluded the state 
equation for suspension feeders in the Chesapeake Bay Eutrophi- 
cation Model, and it combined macrobenthic suspension and 
deposit feeders into a single state equation. The benthic suspension 
feeder equation was not included in our model because its 
construction was based on large bivalve suspension feeders [28] 
that are generally rare in the lower Rappahannock River [43]. 
Based on the community composition of macrobenthos collected 
during the summers of 2007 and 2008 [36], a single governing 
equation for macrobenthos was assumed to be sufficient to model 
macrobenthic response to DO concentrations in the lower 
Rappahannock. Many of the macrobenthos in our samples 
demonstrate both suspension and deposit feeding traits (Table 1), 
and no definable difference in response to hypoxia was observed 
from either group in terms of a change in biomass [36]. The 
combination of these two groups also maintained our goal of 
keeping the model as simple as functionally possible. 

Water quality variables (DIN, POC, and DO) were obtained 
from daily interpolations of Chesapeake Bay Water Quality 
Monitoring Program data from 1985 to 2001 [44], with the 
exception of photosynthetically active radiation (PAR) and water 
temperature. Daily PAR and water temperature were forced using 
equations derived by Wetzel and Neckles [45] for lower 
Chesapeake Bay. 



2.4 Adaptations to the Original Model 

Some specific changes were made to the original governing 
equations of the Chesapeake Bay Eutrophication Model (ad- 
dressed below). The formulation that represented the response of 
zooplankton mortality to hypoxia (Mz) was altered; in our model, 
if DO concentration was less than 2 mg DO 1 1 then: 



Mz = MZERO z * 1 



DOREF \ 
DOCRITz) 



(4) 



in which: 

Mz = hypoxic mortality of zooplankton group Z (d ') 
MZEROz = mortality at zero dissolved oxygen concentration 
(d" 1 ) 



DOREF 



dissolved 



oxygen concentration when DO < 
DOCRIT z , otherwise 2 (mg DO l" 1 ) 

DOCRITz = threshold below which dissolved-oxygen-induced 
mortality occurs, this value equals 2 (mg DO 1 '). 

In the original equation DOCRITz was always 2 mg O z 1 and 
DOREF was the dissolved oxygen concentration when DO < 
DOCRITz, otherwise it was zero. However, during model 
simulation this resulted in a linear increase in modeled zooplank- 
ton population through time, therefore, the formulation was 
amended to equation 4, where the DOREF was the dissolved 
oxygen concentration when DO < DOCRITz, otherwise it was 
2 mg DO V\ 

The parent Chesapeake Bay Eutrophication model simulates 
three fractions of sediment organic carbon, a labile, semi-labile, 
and refractory pool. In the original version of the macrobenthos 

model, the following portion of Eq. 3, a* ( "To*) * 



PLOS ONE | www.plosone.org 



4 



December 2013 | Volume 8 | Issue 12 | e84140 



Modeling the Effect of Hypoxia on Macrobenthos 



POC * k mn \ * B], was computed twice, once for the labile and 
once for the semi-labile carbon pool. Since our model was not 
coupled to a larger eutrophication model, in the interest of 
maintaining simplicity we computed this term in Eq. 3 once using 
total sediment POC from field measurements. 

2.5 Rappahannock Function Relating Biomass to Hypoxia 

In the original Eutrophication model, the impact of DO 
concentration on macrobenthos respiration (r), ingestion (Ig), and 
mortality (m) was represented by the sigmoid equation Z where: 



7 DO g x-DO \ 

i+e \ D °gx-DO qx ) 

where: 

DOj, = DO at which macrobenthos function is 50% of 
maximum. 

T)O qx = DO at which macrobenthos function is 25% of 
maximum. 

In the case of hypoxia-induced mortality, <^is then used in Eq. 5 
above. The sigmoid equation that represents £was not supported 
by any data. Our analysis of Chesapeake Bay field data and results 
from Seitz et al. [24] suggests a sigmoid relationship between DO 
and macrobenthic biomass, as mortality occurs below rather 
sharply defined oxygen concentrations. Therefore, a sigmoid 
function was empirically constrained using macrobenthos data 
collected from the Rappahannock River during the summers of 
2007 and 2008 (Fig. 1), and a function (Z) was fit to the data to 
model the impact of DO concentration on macrobenthicbiomass: 

1 97 

It is important to note that these relationships are quantitative 
representations of the actual various underlying causal factors 
between DO and macrobenthic biomass. In our model, Z 1% still 
used to model the impact of DO concentration on r and Ig (see 
[38]), but Z' is use d to more accurately model the impact of DO 
on macrobenthic biomass, replacing Z m Eq. 5 and providing an 
empirically constrained method for estimating hypoxia-induced 
mortality. Equation 6 was normalized (0 to 1, dimensionless) by 
replacing the numerator with 1, such that Equation 7 was the 
version applied in our model: 

Z = ( DO-1.35 ) 

where: 

DO = dissolved oxygen concentration (mg 1 ). 

2.6 Model Verification and Simulation Analysis 

A single model run encompassed a period of 365 days with a 
time step of one calculation per day. Model output for the three 
governing equations were verified using Chesapeake Bay Benthic, 
Water Quality, and Plankton Monitoring Program data from 
stations LE3.2, LE3.4, and LE3.6 which are located approxi- 
mately 25 km, 19 km, and 2 km from the mouth of the 
Rappahannock River, respectively. While the stations used for 
verification differ spatially, the physical dynamics over this 
relatively small spatial scale of the lower Rappahannock River 



are not significantly different [44]. The macrobenthic output was 
verified with data from stations LE3.2 and LE3.4, and data from 
station LE3.6 were used to verify the phytoplankton and 
zooplankton equations; three stations were used due to data 
availability. Model output was visually compared to observations 
to determine validity [46], and by assessment of mean percent 
error (MPE) and root mean square deviation (RMSD) [32]. 

Predicted phytoplankton biomass was verified against the mean 
annual cycle of chlorophyll-a concentrations computed using 
Chesapeake Bay Program data from 1985 to 2001. To reflect 
mean water column conditions, this annual cycle was based on 
averages of surface and bottom concentrations; differences 
between the surface and bottom were usually less than 10 mg 
m 3 . Predicted zooplankton biomass was verified using computed 
average annual cycles of combined micro- and mesozooplankton 
biomass using Chesapeake Bay Program data from 1985 to 2001. 
Depth-integrated zooplankton counts from the Plankton Moni- 
toring Program were converted to biomass using species-specific 
carbon contents [47]. Predicted macrobenthic biomass was 
verified using Bay Program data from 1992, using site LE 3.4 
for verification under normoxia and site LE3.2 for verification 
under hypoxia. Year 1992 was chosen at random from years 
1985-2001. The approach of using a randomly selected year for 
macrobenthos verification instead of averaged annual cycles was 
due to the stochastic nature of benthic data collection through 
time; data were collected at variable dates each year by the 
Chesapeake Bay Program. Following verification, sensitivity 
analyses were conducted for the phytoplankton, zooplankton, 
and macrobenthic state variables by adjusting selected parameters 
that directly impacted growth or loss (i.e. consumption or 
predation parameters). Maximum photosynthetic rate (P^rri), 
phytoplankton settling velocity (Wd), and predation rate on algae 
(Phtl) were tested for the phytoplankton state variable; predator 
biomass and clearance rate (PHTlz) for the zooplankton state 
variable; and assimilation efficiency for carbon (a) and ingestion 
limitation [Kj) for the macrobenthic state variable. All parameters 
tested in sensitivity analysis were adjusted at an increment of 
±20% and the relative percent difference from the standard run 
was calculated for each. Parameters with percentage differences 
greater than 10% were deemed to be sensitive parameters [48]. 

A set of simulations analyses were then conducted, adjusting 
DO concentration to model the effect of the severity and duration 
of hypoxia on the three modeled state equations (Table 2). DO 
concentrations were gradually adjusted during a period of 3 days 
to avoid artificial rapid changes in the model. Simulations SI to S3 
focused on the sustained duration of hypoxia. S4 simulated the 
development of hypoxia during the neap-spring tidal cycle with 
DO cycling for 14 days. Verified model output from the state 
equations (i.e. the base model results during normoxia after model 
verification and sensitivity analyses; Fig. 2, 3, and 4) was used as a 
baseline for comparison against simulations SI to S4. This was 
conduct through visual observation and comparison of the means 
and standard deviations between the verified model output and 
the simulations (e.g. phytoplankton, macrobenthos, and zooplank- 
ton output in SI was compared to the models in Fig. 2, 3, and 4, 
respectively). Simulations S5 to S9 modeled the severity of hypoxia 
from 0.0 to 2.0 mg 0 2 1 _1 with a 0.5 mg O z l -1 step and constant 
duration of 60 days. Comparison of means and standard 
deviations were used to justify "significant" differences. DO 
concentration was forced for each model simulation, with 
consideration of historical data of DO dynamics in the 
Rappahannock River [44] . 
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0.0 

0 73 146 219 292 365 
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Figure 2. Verification of the phytoplankton state variable. Left y-axis: the gray solid line represents modeled phytoplankton biomass and the 
black dots denote mean observed phytoplankton biomass from site LE3.6 of the Chesapeake Bay Water Quality Monitoring Program from 1 985-2001 . 
Right y-axis: the grey dashed line indicates the dissolved oxygen concentration. Modeled phytoplankton biomass matched the trend and magnitude 
of observed phytoplankton biomass in 1992. 
doi:1 0.1 371 /journal.pone.00841 40.g002 



Model 




146 219 

Julian Day 



Figure 3. Verification of the macrobenthos state variable during normoxia. Left y-axis: the gray solid line represents modeled 
macrobenthos biomass and the black dots denote observed macrobenthos biomass from site LE3.4 of the Chesapeake Bay Benthic Monitoring 
Program in 1992. Right y-axis: the grey dashed line indicates the dissolved oxygen concentration. Modeled macrobenthos biomass matched the 
trends and magnitude of observed macrobenthos biomass in 1992. 
doi:1 0.1 371 /journal.pone.00841 40.g003 
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Figure 4. Verification of the zooplankton state variable. Left y-axis: the gray solid line represents modeled zooplankton biomass and the black 
dots denote mean observed zooplankton biomass from site LE3.6 of the Chesapeake Bay Water Quality Monitoring Program from 1985-2001. Right y- 
axis: the grey dashed line indicates the dissolved oxygen concentration. Modeled zooplankton biomass approximated the observed annual cycle and 
matched the magnitude of observed zooplankton biomass in 1992. 
doi:1 0.1 371 /journal. pone.00841 40.g004 



Results 

Modeled phytoplankton biomass (Pj matched the general trends 
in observed biomass in the lower Rappahannock River, with 
blooms and declines consistent with the data from 1992 (Fig. 2). 
The model also reproduced the approximate magnitude of 
phytoplankton biomass with MPE of 13.3% and RMSD of 
0.18 g C m . Modeled macrobenthic biomass (B) also matched 
the patterns in observed biomass in 1992 at both a normoxic site 
(Fig. 3) and a hypoxic site (Fig. 5). The model accurately portrayed 
the temporal dynamics in macrobenthic biomass with respect to 
DO concentration, as well as the magnitude of biomass with MPE 
of 12.1 % and RMSD of 0.56 during normoxia and MPE of 15.7% 



and RMSD of 0.22 g C m 2 during hypoxia. Modeled 
zooplankton biomass matched the magnitude of observed biomass 
in 1992, but did not follow the trends in the observed data (Fig. 4), 
with MPE of 4.5% and RMSD of 0.38 g C m" 2 . Modeled 
zooplankton dynamics appeared to be delayed when compared to 
observed increases and decreases in zooplankton biomass. This 
discrepancy may be due to the combination of micro- and 
mesozooplankton into a single state variable. Additionally, there 
were noted calibration difficulties for zooplankton in the parent 
eutrophication model [49], which may have extended into our 
adaptation of that model. 

Discrepancies between modeled and observed zooplankton 
biomass could be expected to have detrimental effects on modeled 



Table 2. Model simulations. 





Simulations Hypoxia Duration (d) 


Julian Day 


Ordinal Date 


DO (mg 0 2 I 1 ) 


S1 120 


148-267 


May 28- Sept 24 


0.5 


S2 60 


178-237 


Jun 26- Aug 26 


0.5 


S3 30 


191-222 


Jul 11- Aug 11 


0.5 


S4 14 d intervals 


162-176, 190-204, 218- 
232, 246-260 


Jun 1 1-25, Jul 9-23, Aug 6-20, 
Sept 3-17 


0.5 


S5 60 


178-237 


Jun 26- Aug 26 


2.0 


S6 60 


178-237 


Jun 26- Aug 26 


1.5 


S7 60 


178-237 


Jun 26- Aug 26 


1.0 


S8 60 


178-237 


Jun 26- Aug 26 


0.5 


S9 60 


178-237 


Jun 26- Aug 26 


0.0 


Ecosystem model simulations with varying hypoxic duration and severity. SI to S4 (Figures 6 and 7) modeled hypoxic duration at a constant concentration of 0.5 mg 0 2 
r 1 ; S4 modeled intermittent hypoxia (hypoxia occurring every 14 days on a neap/spring tidal cycle); S5 to S9 modeled the affect of hypoxic severity at a constant 



duration of 60 days. 
doi:1 0.1 371 /journal.pone.0084140.t002 
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0.00 J — — i — — I- 0 
0 73 146 219 292 365 

Julian Day 

Figure 5. Verification of the macrobenthos state variable during hypoxia. Left y-axis: the gray solid line represents modeled macrobenthos 
biomass and the black dots denote observed macrobenthos biomass from site LE3.2 of the Chesapeake Bay Benthic Monitoring Program in 1992. 
Right y-axis: the grey dashed line indicates the dissolved oxygen concentration. Modeled macrobenthos biomass matched the trends and magnitude 
of observed macrobenthos biomass in 1992. 
doi:1 0.1 371 /journal.pone.00841 40.g005 



phytoplankton biomass, given the significant influence of zoo- 
plankton predation rate on algae (see Table 3 below). However, 
modeled phytoplankton biomass successfully reproduced the 
observations (Fig. 2). Additionally, phytoplankton biomass in the 
lower Rappahannock River is 3.3-14.3 times greater (mean = 6.5) 
than zooplankton biomass (Fig. 2 and 4), which indicates that 
phytoplankton will have the dominant effect on modeled 
macrobenthic biomass. The successful simulation of phytoplank- 
ton biomass, simulation of zooplankton biomass in the correct 
range (Fig. 4), the overall much greater biomass of phytoplankton 
compared to zooplankton, and the accurate simulation of 
macrobenthic biomass under both normoxia (Fig. 3) and hypoxia 
(Fig. 5), confirm the ability of the model to predict the dynamics of 
macrobenthic biomass in response to hypoxia in further simulation 
analyses. 

Sensitivity analyses were conducted on model parameters for 
each state variable (Table 3). The model was sensitive to a majority 

Table 3. Sensitivity analysis. 



of tested parameters, with the phytoplankton state variable 
sensitive to all tested parameters. The zooplankton state variable 
was found to be sensitive to an increase and decrease in predator 
biomass and clearance rate. The macrobenthic state variable was 
sensitive to assimilation efficiency for carbon and insensitive to 
ingestion limitation. The sensitivity of these parameters is a 
problem in the parent model that we did not attempt to address. 
The goal of this manuscript was not to resolve nuances in the 
original model, but to take a simple approach and constrain the 
effect of DO concentration on macrobenthic biomass, within the 
context of the more complex Chesapeake Bay Eutrophication 
Model. 

Simulations were run assessing the impact of hypoxic duration 
on phytoplankton, zooplankton, and macrobenthic biomass. 
Macrobenthic biomass began to decrease as scenario SI, which 
simulated a hypoxic duration of 120 days, approached hypoxia 
(Fig. 6A). At the start of hypoxia in SI the steady decrease in 





State variable 


Parameter 


-20% 


+20% 


Average RMS 


% Diff -20% 


% Diff +20% 


Phytoplankton 


P"m 


0.032 


0.041 


0.037 


12.7* 


11.8* 




Wa 


0.041 


0.032 


0.037 


11.5* 


12.2* 




Phtl 


0.044 


0.030 


0.037 


18.8* 


18.6* 


Zooplankton 


PHTIz 


0.001 


0.003 


0.002 


11.8* 


13.1* 


Macrobenthos 


a. 


0.118 


0.369 


0.192 


38.8* 


91.8* 




K, 


0.191 


0.194 


0.192 


0.8 


0.5 



Results of sensitivity analysis for phytoplankton, zooplankton, and macrobenthic state variables. The root mean square (RMS) was the average variance in the base state 

equation over all time steps of a single year (n = 365). RMS values are shown for ±20% variation for each state variable by parameter. The model was deemed to be 

sensitive when % difference exceeded 1 0% (Ragone-Calvo et al. 2001 ). Parameters: Maximum photosynthetic rate {P b m), phytoplankton settling velocity [Wa), predation 

rate on algae {Phtl), predator biomass and clearance rate {PHTIz), assimilation efficiency for carbon (a), and ingestion limitation {Kj). 

'Denotes model sensitivity. 

doi:1 0.1 371 /journal.pone.0084140.t003 
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0 73 146 219 292 365 

Julian Day 

Figure 6. Simulations of hypoxic duration and macrobenthic 
response. Simulated macrobenthos biomass (8) under hypoxic 
durations of (A) 120, (B) 60, and (C) 30 days; SI, S2, and S3, respectively. 
Left y-axis: the gray solid line represents modeled macrobenthos 
biomass and the black dots denote observed macrobenthos biomass 
from site LE3.4 of the Chesapeake Bay Benthic Monitoring Program in 
1992. Right y-axis: the grey dashed line indicates the dissolved oxygen 
concentration. S was influenced by the duration of hypoxia. 
doi:1 0.1 371 /journal.pone.00841 40.g006 

macrobenthic biomass accelerated to an immediate collapse in 
biomass that lasted the duration of the hypoxic event, with a 
temporal trend in macrobenthic biomass different from that 
modeled under normoxic conditions. Macrobenthic biomass 
began to respond and increase before DO concentrations in SI 
became normoxic. A few days after hypoxia ended, macrobenthic 
biomass had increased to above pre-hypoxia levels. Similar trends 
were observed in S2 and S3 involving hypoxic durations of 60 and 
30 days (Fig. 6B and C), with the length of the collapse the main 
difference in macrobenthic biomass. Simulations with shorter 
durations of hypoxia resulted in less time with macrobenthic 
biomass near 0 g C m -2 . Macrobenthic biomass of S2, with a 
mean of 4. 1 g G m (SD = 3.1), and S3, with a mean of 4.4 g C 
m -2 (SD = 2.8), were different than the macrobenthic biomass 
model under normoxia, with a mean of 3.5 g G m" 2 (SD= 1.5). 

In simulation S4, with intermittent hypoxia, macrobenthic 
biomass decreased at the onset of hypoxia and remained near zero 
for the duration of hypoxia (Fig. 7). During the 14-day intervals 
when DO concentration was normoxic, macrobenthic biomass 
began to increase but decreased back near zero with the onset of 
hypoxia. Macrobenthic biomass fluctuated through this pattern 
throughout the hypoxic simulation. A few days after hypoxia 
abated permanently, macrobenthic biomass began to increase to 
biomass levels greater than those observed pre-hypoxia. 

Modeled phyto- and zooplankton biomass responded oppositely 
during hypoxic simulations (Fig. 8). As DO concentrations began 
to decline, phytoplankton biomass initially decreased, however, 
with the onset of hypoxia, phytoplankton biomass increased. The 
length of increased phytoplankton biomass was dependent on the 



duration of hypoxia, with a longer duration of hypoxia resulting in 
higher overall phytoplankton biomass, and to some extent an even 
greater magnitude of phytoplankton biomass. Hypoxia had the 
reverse effect on zooplankton biomass. As DO concentration 
decreased to hypoxic levels, zooplankton biomass initially 
increased and then declined to near 0 g C m ' 5 . The length of 
time that zooplankton biomass stayed near 0 g C m 3 was 
dependent on the duration of hypoxia; lengthy durations of 
hypoxia coincided with longer durations of reduced zooplankton 
biomass. Modeled zooplankton biomass did not respond to 
increased DO concentration until days to weeks after hypoxia 
ended. 

The effect of hypoxic severity on macrobenthic biomass was 
tested by adjusting DO concentration between 2.0 and 0.0 mg 0 2 
l -1 in increments of 0.5 mg O z 1 _1 . Macrobenthic biomass was 
not substantially different between simulations when compared 
during the full year (Table 4). The mean biomass of these four 
simulations was approximately equal, with large standard devia- 
tions. However, these simulations had similar macrobenthic 
biomass during normoxia, as no parameters were changed; 
macrobenthic biomass did not differ between simulations until 
DO concentrations became hypoxic. Therefore, simulations S5 to 
S9 were analyzed starting at the onset of hypoxia on day 178 
through the end of hypoxia on day 237 (Table 5). S5 had 
appreciably higher macrobenthic biomass than S6 to S9. Mean 
biomass in S5, which depicted 60 days of hypoxia at 2.0 mg 0 2 
1 , was >5 times higher than simulations that modeled hypoxia 
at DO concentrations of 0.0— 1.0 mg 0 2 1 , and 3 times higher 
than S6. In S6 when the DO was 1.5 mg 0 2 1 1 macrobenthic 
biomass was lower than S5 and >2 times higher than S7 to S9. 
There were negligible differences in biomass at DO concentrations 
of 1.0 mg 0 2 1 1 and lower. The mean biomass in S7, S8, and S9 
was approximately similar with large standard deviations. 

Discussion 

The sigmoid function (£) applied to this ecosystem model 
resulted in clear changes in simulated macrobenthic biomass (B) 
during various hypoxic scenarios. In our model, hypoxic duration 
resulted in prolonged reductions of macrobenthic biomass relative 
to the length of hypoxia, with the model suggesting near 
defaunation (macrobenthic biomass equal to 0 g C m~ 2 ) during 
the 120, 60, and 30 day hypoxic scenarios (simulations 1, 2, and 3, 
Fig. 6) at a DO concentration of 0.5 mg 0 2 1 _1 . This is in good 
agreement with laboratory and field studies that show duration 
and severity of hypoxia to impact benthic communities (see 
[11,26]). Periods of prolonged hypoxia have been observed 
previously in Chesapeake Bay and elsewhere [8,50]. The deep 
trough of the mainstem Bay experiences sustained seasonal 
hypoxia year after year [15]. While some species have documented 
resistance to hypoxia [7], during extended periods of hypoxic 
exposure (~40 days) even the most tolerant of species experience 
total mortality [51]. 

After DO levels in our model returned to normoxia, macro- 
benthic biomass recovered greater than pre-hypoxic levels. This 
was unexpected as conditions in Chesapeake Bay that stimulate 
ecological production in the pre-hypoxic spring differ in the post- 
hypoxic fall [8] . Increased nutrient run-off from the spring freshet 
promotes plankton production; the particulate organic matter 
from these blooms eventually settles to the bottom promoting 
benthic growth [52]. Large plankton blooms seen in the spring are 
noticeably absent in the fall, and with less primary production one 
would expect the rate of recovery of macrobenthic biomass to be 
less in the fall than in the spring. However, data from the 
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Figure 7. Simulations of intermittent hypoxic duration and macrobenthic response. Simulated macrobenthos biomass (B) under 
intermittent hypoxia (14 days hypoxic and 14 days normoxic) representing effect of neap/spring tide cycle. Left y-axis: the gray solid line represents 
modeled macrobenthos biomass and the black dots denote observed macrobenthos biomass from site LE3.4 of the Chesapeake Bay Benthic 
Monitoring Program in 1992. Right y-axis: the grey dashed line indicates the dissolved oxygen concentration on the right y-axis. 8 responded to 
modeled intermittent hypoxia. 
doi:1 0.1 371 /journal.pone.00841 40.g007 




Julian Day 

Figure 8. Simulations of hypoxic duration and phyto- and zooplankton response. Simulated phytoplankton biomass under hypoxic 
durations of (A) 1 20, (B) 60, and (C) 30 day, and zooplankton biomass under hypoxic durations of (D) 1 20, (E) 60, and (F) 30 days. Left y-axis: the gray 
solid line represents modeled biomass and the black dots denote mean observed biomass from site LE3.6 of the Chesapeake Bay Water Quality 
Monitoring Program from 1985-2001. Right y-axis: the grey dashed line indicates the dissolved oxygen concentration on the right y-axis. Modeled 
phytoplankton and zooplankton biomass responded inversely to hypoxic simulations. 
doi:1 0.1 371 /journal.pone.00841 40.g008 
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Table 4. Response of macrobenthic biomass to hypoxic 
severity (365 d). 





Simulation DO (mg 0 2 1 ') 


Julian Day 


Mean Biomass (g) 


S5 2.0 


1-365 


3.20 (1.84) 


S6 1.5 


1-365 


3.04 (2.07) 


S7 1.0 


1-365 


2.88 (2.07) 


S8 0.5 


1-365 


2.88 (2.07) 


S9 0.0 


1-365 


2.88 (2.07) 


Comparison of macrobenthos biomass to hypoxic severity over a full year. 



Macrobenthos biomass was not measurably different between simulations over 
a full year. 

doi:1 0.1 371/joumal.pone.0084140.t004 



continuously monitored sites used to derive Z' indicated macro- 
benthic production can increase to pre-hypoxic levels a few weeks 
post hypoxia, suggesting the macrobenthic biomass increases we 
observed post-hypoxia may be reasonable [36]. A decrease in 
predation pressure may also have contributed to the observed 
increase in macrobenthic biomass. In an ecological context, 
decreases in temperature during the fall and winter months post- 
hypoxia lower metabolic demand, and thus the need of predators 
to obtain food [53]. In our model this was accounted for through 
Arrhenius temperature dependencies in the predation formulation, 
which is a simple but accurate formula for the temperature 
dependence of reaction rates [54]. This however, suggests that 
during summer hypoxia, the loss of macrobenthic biomass and its 
impacts to higher consumers such as epibenthic predators and 
demersal fish is not compensated for by the increase in 
macrobenthic biomass post-hypoxia due to a change in metabolic 
demand and activity. 

Modeled intermittent hypoxia (S4, Fig. 7) resulted in reduced 
macrobenthic biomass during hypoxic periods but some recovery 
during normoxia. This cycling continued throughout the inter- 
mittent series of hypoxic events. Given the severity at which the 
DO concentration was set for this simulation (0.5 mg O2 1~ ), it is 
not surprising that macrobenthic biomass decreased to the 
observed level. Neubauer [55] found in situ decreases in 
instantaneous macrobenthic production coincided with a hypoxic 
event, however observed no cyclical pattern between macro- 
benthic production and hypoxia. Mean DO concentration during 
Neubauer's [55] hypoxic event were —2.1-2.5 mg O z 1 , which 
would be considered mild hypoxia for macrobenthos [56]. 
Further, Neubauer [55] experienced large recruitment events for 
some macrobenthic species during hypoxia at the time of his study. 
In contrast, a cyclical pattern between hypoxia and macrobenthic 
production was observed in data from the continuously monitored 
sites used to derive Z'- DO concentrations at our in situ sites were 
as low as ~0.02 mg 0 2 1 1 on two separate occasions, and the 
similarities in hypoxic severity between our model and our in situ 
observations may account for analogous results [36]. 

As with simulations SI, S2, and S3, the growth of macrobenthic 
biomass during simulation S4 was very rapid with the return of 
normoxia, an observation that was also documented in situ [36] . A 
number of the dominant macrobenthic species in the lower 
Rappahannock spawn and recruit throughout the summer [57- 
63] . The larvae of these species are found distributed throughout 
the water column, and this pattern does not change in response to 
low DO [64]. When favorable conditions return, the planktonic 
larvae of these benthic species are available to settie in an organic- 
rich environment with few competitors, which may explain why 



Table 5. Response of macrobenthic biomass to hypoxic 
severity (60 days). 




Simulation 


DO (mg 0 2 1 1 


Julian Day 


Mean Biomass (g) 


S5 


2.0 


178-237 


4.07 (0.46) 


S6 


1.5 


178-237 


1 .77 (0.69) 


S7 


1.0 


1 78-237 


0.76 (1.15) 


S8 


0.5 


178-237 


0.62 (1.15) 


S9 


0.0 


1 78-237 


0.62 (1.15) 



Comparison of macrobenthos biomass to hypoxic severity over a partial year, 
covering the 60 day time-frame of simulated hypoxia. Macrobenthos biomass 
was measurably variable during the 60 day assessment. Mean biomass is shown 
with ±1 SD in parentheses. 
doi:1 0.1 371/joumal.pone.0084140.t005 



quick benthic recruitment post-hypoxia is observed. Modeled 
macrobenthic biomass also began recovering ~2-3 days before 
hypoxia abated, which represents the time frame that DO was 
increasing from 0.5 to 2.0 mg O z 1 , indicating a sensitivity to 
hypoxic severity. Macrobenthos have been previously documented 
to recruit before the DO concentration is above 2.0 mg 0 2 1 _ , 
although this low level of DO does delay settlement of some 
hypoxia sensitive taxa [64]. In our intermittent simulation, once 
macrobenthic biomass reached a level where it appeared 
sustainable, hypoxia returned and biomass was reduced to near 
zero again. 

Hypoxic severity had a measureable impact on macrobenthic 
biomass with higher biomass in simulations with less severe 
hypoxia. Studies have shown the severity of hypoxia to affect the 
response of benthic communities; the more severe the hypoxia, the 
greater the impact on the benthos, directly and indirecdy [7]. 
Directly, benthic species vary in their tolerances to low DO 
concentrations [11], and as the severity of hypoxia increases 
towards anoxia, sensitive species die-off, usually due to asphyxi- 
ation, decreasing the diversity of the affected area and overall 
biomass [7] . Additionally, the toxic compound hydrogen sulfide is 
present in severely hypoxic sediments, and has been suggested as a 
mechanism contributing to macrobenthic mortality during hyp- 
oxia [65]. Indirectly, DO concentrations can positively and 
negatively affect benthic predation; Nestlerode and Diaz [66] 
showed that benthos may actually have a refuge from predation 
under mild hypoxic conditions, and Brante and Hughes [67] 
demonstrated that hypoxia reduced Carcinus maenas predation on 
mussels. If macrobenthos are able to avoid mortality via 
asphyxiation, such actions during prolonged hypoxic events could 
also indirectly lead to starvation [7] . 

During model simulations, there was no increase in macro- 
benthic biomass as DO concentrations declined towards hypoxia. 
In our model B: 



/> = /?' 



DO 

DO + K DO 



(8) 



in which:/? = the predation rate before considering hypoxic 
tStc\sK D0 - predation DO half-saturation. 

accounts for the predation rate on macrobenthos and denotes 
predation rate as a function of temperature and DO concentra- 
tion, where predation increases with increasing DO. During mild 
hypoxia predators may not effectively prey upon benthos, and 
hypoxia tolerant benthos would survive and maintain their 
biomass. The model excludes the effect of increasing predation 
during hypoxia, an artifact of the parent model. Seitz et al. [68] 
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and Long and Seitz [21] showed that epibenthic predators and 
demersal fish can at times capitalize on stressed benthos during 
mild hypoxic events. As oxygen concentrations become lethal, 
stressed macrobenthos extend their appendages and bodies out of 
the sediment in an attempt to escape severely hypoxic conditions 
below the sediment-water interface [56]. Opportunistic mobile 
predators have been shown to re-enter hypoxic areas and prey on 
exposed macrobenthos during mild hypoxia [69,70]. 

Zooplankton biomass (M) was negatively impacted by hypoxia 
directly, causing zooplankton biomass to be drastically reduced. 
Marcus et al. [71] considered the effect of reduced DO 
concentration on the survival and population dynamics of 
zooplankton, demonstrating the deleterious effect hypoxia has on 
zooplankton population and community dynamics. In our model, 
phytoplankton biomass was sensitive to predation by zooplankton, 
with a 20% increase and 20% decrease in algal predation resulting 
in an 18% difference in phytoplankton biomass (Table 3). As a 
result, phytoplankton biomass (P) in our model was indirecdy 
positively influenced by the onset of hypoxia, due to the release of 
phytoplankton from zooplankton grazing pressure. 

Conclusions 

Macrobenthos data from the lower Rappahannock River were 
used to derive a sigmoid relationship, to model the effect of DO 
concentration on macrobenthic biomass (B). Z' was then used to 
empirically constrain a formulation for hypoxia-induced mortality 
in a biomass-based ecosystem model and used to assess the impact 
of hypoxia on macrobenthic biomass, while including the 
important biological interactions that occur through benthic- 
pelagic coupling. Z' is a useful tool in that it can be applied to 
existing models to simulate the impact of hypoxia on the 
macrobenthos, and the methods used to derive Z' can be applied 
to other systems to develop site specific and species specific 
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parameterizations of Z'- Further, in terms of application, the 
simplicity of our model makes it easy to implement, interpret, and 
update, while at the same time the approach developed here can 
be readily incorporated into more complex models. 

From our modeling efforts we found that the duration and 
severity of hypoxia negatively affected macrobenthic biomass; 
longer durations and greater hypoxic severity resulted in less 
biomass. Further, our model suggests that post hypoxia, macro- 
benthic biomass can return to pre-hypoxic levels, which implies a 
level of resiliency in the macrobenthic community to hypoxia. 
However, it is important to remember that the loss of macro- 
benthic biomass occurs at a critical time when energy demands of 
epibenthos and demersal fish, predators of macrobenthic organ- 
isms, are at their highest. It is unknown if the rebound in 
macrobenthic biomass observed in our model post-hypoxia would 
be enough to compensate for the loss in biomass during hypoxia. 

The ecological importance of macrobenthos to estuarine 
systems underlies the significance in understanding processes that 
positively and negatively impact this group. The ability to 
accurately model the influence of low DO and strength of 
interactions between ecosystem components will improve our 
understanding of the impacts of hypoxia and provide a holistic 
view of a major anthropogenic stressor on ecosystem functioning. 
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